In [33]:
# Import libraries
import pandas as pd
import os
import numpy as np
import matplotlib.pyplot as plt
import altair as alt
import altair_saver 
import warnings
import geopandas as gpd
import cartopy.io.img_tiles as cimgt
import cartopy.feature as cfeature
import xarray as xr
from shapely.geometry import Point
import folium
from folium.plugins import MarkerCluster
In [34]:
os.chdir('C://Users/pmarshal/Documents/Climate-Summary')
os.getcwd()
Out[34]:
'C:\\Users\\pmarshal\\Documents\\Climate-Summary'

Precipitation¶

The charts below show precipitation trends in the watersheds. You will see plots with annual, monthly, and daily precipitation amounts. Atmospheric river events are also highlighted. The reference station for these plots is the Lower Capilano Fire Weather station, on the east side of the Capilano Reservoir. Precipitation amounts tend to be greater further to the north.

In [35]:
precip = pd.read_csv('data/precipitation_2023.csv')
#df['year'] = pd.to_datetime(df['year'], format='%Y')
#df['date'] = pd.to_datetime(df['date'], format='%Y-%m-%d')
#precip = precip.drop(labels=[25,26], axis=0)
#precip.head()
In [36]:
precip2 = pd.read_csv('data/precipitation_2024.csv')
#df['year'] = pd.to_datetime(df['year'], format='%Y')
#df['date'] = pd.to_datetime(df['date'], format='%Y-%m-%d')
#precip = precip.drop(labels=[25,26], axis=0)
#precip2.head()
In [37]:
precip2 = precip2.sort_values("precip")
#precip2.head()
In [38]:
title_pc = alt.TitleParams(
   text='2024 Annual Precipitation - Metro Vancouver',
#   subtitle="Gray area is 5th to 95th percentile. Red dashed line is 20-year mean",
   anchor='middle',
   fontSize=14,
   fontWeight='bold')

bar = alt.Chart(precip2, title=title_pc).mark_bar().encode(
    x=alt.X('Site:O', title=None, sort= 'y'),
    y=alt.Y('precip:Q', title='Total Annual Precipitation (mm)'),
    color=alt.Color('precip:Q', scale=alt.Scale(scheme='viridis', reverse=True), legend=None)
)


text2 = bar.mark_text(
    align='left',
    baseline='middle',
    dx=5,  # Nudges text to right so it doesn't appear on top of the bar
    angle=270
).encode(
    text='precip:Q'
)
(bar + text2).properties(width=600)
Out[38]:
In [39]:
# Specify the CRS (WGS84, EPSG:4326)
crs = "EPSG:4326"

# Create the geometry with specified CRS
geometry = [Point(xy) for xy in zip(precip.Long, precip.Lat)]
precip_points = gpd.GeoDataFrame(precip, geometry=geometry, crs=crs)

#precip_points.head(10)

The map below shows the total annual precipitation for a number of weather stations in the region. You can click on the stations to see the total precipitation for that site.

In [40]:
desired_min = 1300
desired_max = 5200
# Create a base map
#map_center = [precip['Lat'].mean(), precip['Long'].mean()]
mymap = folium.Map(location=[49.426, -123.000], zoom_start=9)
#mymap = folium.Map(location=map_center, zoom_start=10)

# Add markers with different colors based on precipitation values
marker_cluster = MarkerCluster().add_to(mymap)

# Define a color scale for precipitation values
color_scale = folium.LinearColormap(colors=['#FDE725', '#6DCC73', '#26828E', '#3E4A89', '#440154'], vmin=desired_min, vmax=desired_max)

# Add markers with different colors based on precipitation values
for index, row in precip.iterrows():
    color = color_scale(row['precip'])
    popup_text = f"Site: {row['Site']}<br>Precipitation: {row['precip']} mm"
    folium.CircleMarker(location=[row['Lat'], row['Long']], radius = 8, color=color, fill=True, fill_color=color, fill_opacity=0.7, popup=popup_text).add_to(mymap)

# Add the color scale to the map
color_scale.add_to(mymap)

# Save the map to an HTML file
mymap
Out[40]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [41]:
# Import dataset
qd_pc = pd.read_csv('data/QD1_temps.csv', parse_dates=['date']) 
qd_pc['wtr_day_yr'] = pd.to_datetime('2023-01-01') + pd.to_timedelta(qd_pc['wtr_yr_day'] - 93, unit='D')
qd_pc = qd_pc.dropna()
qd_pc['year'] = pd.DatetimeIndex(qd_pc['date']).year
qd_pc['month'] = pd.DatetimeIndex(qd_pc['date']).month
qd_pc['day'] = pd.DatetimeIndex(qd_pc['date']).day
qd_pc['DOY'] = pd.DatetimeIndex(qd_pc['date']).dayofyear
qd_pc['quarter'] = pd.DatetimeIndex(qd_pc['date']).quarter
#qd_pc.head()
In [42]:
#Caluclate statistics
qd_pc['mean'] = qd_pc.groupby('wtr_yr_day')['rain_cum_wyr'].transform('mean').round(2)
qd_pc['max'] = qd_pc.groupby('wtr_yr_day')['rain_cum_wyr'].transform('max').round(2)
qd_pc['min'] = qd_pc.groupby('wtr_yr_day')['rain_cum_wyr'].transform('min').round(2)
qd_pc['std'] = qd_pc.groupby('wtr_yr_day')['rain_cum_wyr'].transform('std').round(2)
qd_pc['sem'] = qd_pc.groupby('wtr_yr_day')['rain_cum_wyr'].transform('sem').round(2)
qd_pc['ci95_hi'] = qd_pc['mean'] + 1.96* qd_pc['sem']
qd_pc['ci95_lo'] = qd_pc['mean'] - 1.96* qd_pc['sem']
qd_pc['ci95_hi'] = qd_pc['ci95_hi'].round(2)
qd_pc['ci95_lo'] = qd_pc['ci95_lo'].round(2)
qd_pc['percentile_5th'] = qd_pc.groupby('wtr_yr_day')['rain_cum_wyr'].transform(lambda x: np.percentile(x, 5))
qd_pc['percentile_95th'] = qd_pc.groupby('wtr_yr_day')['rain_cum_wyr'].transform(lambda x: np.percentile(x, 95))
qd_pc['percent_5th'] = qd_pc.groupby('DOY')['rain_cum_yr'].transform(lambda x: np.percentile(x, 5))
qd_pc['percent_95th'] = qd_pc.groupby('DOY')['rain_cum_yr'].transform(lambda x: np.percentile(x, 95))
#qd_pc.tail()
In [43]:
qd_pc = qd_pc.drop(['mean_TA', 'min_TA', 'max_TA'], axis=1)
In [44]:
last_few = qd_pc.loc[qd_pc['wtr_yr'] >= 2022]
#last_few.tail()
# Export DataFrame to CSV
#last_few.to_csv('output.csv', index=False)
In [45]:
qd_monthly = qd_pc.groupby(['year', 'month'])['rain'].sum().reset_index()
qd_monthly['date'] = pd.to_datetime(qd_monthly[['year', 'month']].assign(day=1))
qd_monthly['mean_rain'] = qd_monthly.groupby('month')['rain'].transform('mean').round(1)
qd_monthly['percent'] = qd_monthly['rain']/qd_monthly['mean_rain'].round(1)
qd_monthly['quarter'] = pd.DatetimeIndex(qd_monthly['date']).quarter
#qd_monthly['mean'] = qd_pc.groupby(['year', 'month'])['rain'].mean().reset_index()
#qd_monthly.tail(12)
In [46]:
# Define a function to map months to meteorological seasons
def get_season(month):
    if month in [12, 1, 2]:
        return 'Winter'
    elif month in [3, 4, 5]:
        return 'Spring'
    elif month in [6, 7, 8]:
        return 'Summer'
    else:
        return 'Fall'

# Apply the function to create a new 'season' column
qd_monthly['season'] = qd_monthly['month'].apply(get_season)


monthly_2024 = qd_monthly.tail(13)
#monthly_2024
In [47]:
# Now you can group by 'season' and aggregate as needed
seasonal_totals = monthly_2024.groupby('season')['rain'].sum().reset_index()
#seasonal_totals
In [48]:
quarter_sum = qd_monthly.groupby(['year', 'quarter'])['rain'].sum().reset_index()
year_sum = qd_monthly.groupby(['year'])['rain'].sum().reset_index()
quarter_two = quarter_sum.loc[quarter_sum['quarter'] == 3]
#year_sum
In [49]:
rain = qd_pc.loc[qd_pc['year'] == 2024]
rain_mo = rain.loc[rain['month'] == 9]
rain.to_csv('qd_rain.csv', index=False)
#rain.head()

The chart below shows daily precipitation amounts at the Lower Capilano Fire Weather Station. On October 19, 216 mm of rain fell during an atmospheric river event.

In [50]:
day_bar = alt.Chart(rain).mark_bar(color='darkgreen', size=3).encode(
    x=alt.X('monthdate(date):T', title=None),
    y=alt.Y('rain:Q', title='Precipitation (mm)'),
    tooltip=[alt.Tooltip('date', title="Date"), alt.Tooltip('rain', title="Precip (mm)")]
).properties(width=700, title='2024 Daily Precipitation - Lower Capilano Watershed')

day_bar
Out[50]:
In [51]:
# Melt the dataframe to have a 'Days' column indicating whether it's 'Dry' or 'Wet'
#df_melted = pd.melt(merged_df, id_vars=['Year'], var_name='Day Type', value_name='Days')
In [52]:
ar_df = pd.read_csv('data/AR_2023.csv')
ar_df['date'] = pd.to_datetime(ar_df['date'])
ar_df['duration'] = ar_df['duration'].astype(float)
ar_df['precip_cum'] = ar_df['precip_cum'].astype(float)
ar_df['year'] = pd.DatetimeIndex(ar_df['date']).year
#ar_df.tail()

The next plot below breaks down precipitation amounts by month. The box and whiskers plots show the distribution of precipitation each month, the blue or orange bars show this years values, and the plot on the right shows the percent of normal precipitation.

In [53]:
rain_sum = rain.groupby(['quarter'])['rain'].sum().reset_index()
#rain_sum
In [54]:
# Add text annotation
text1 = alt.Chart().mark_text(
    text='Winter = 769mm',  # Replace with your desired text
    align='left',
    size=12,
    baseline='middle',
    dx=0  # adjust the position of the text
).encode(
    x=alt.value(10),  # Adjust the x-position of the text box
    y=alt.value(10)   # Adjust the y-position of the text box
)

text2 = alt.Chart().mark_text(
    text='Spring = 373mm',  # Replace with your desired text
    align='left',
    size=12,
    baseline='middle',
    dx=160  # adjust the position of the text
).encode(
    x=alt.value(10),  # Adjust the x-position of the text box
    y=alt.value(10)   # Adjust the y-position of the text box
)

text3 = alt.Chart().mark_text(
    text='Summer = 197mm',  # Replace with your desired text
    align='left',
    size=12,
    baseline='middle',
    dx=310  # adjust the position of the text
).encode(
    x=alt.value(10),  # Adjust the x-position of the text box
    y=alt.value(10)   # Adjust the y-position of the text box
)

text4 = alt.Chart().mark_text(
    text='Fall = 1103mm',  # Replace with your desired text
    align='left',
    size=12,
    baseline='middle',
    dx=470  # adjust the position of the text
).encode(
    x=alt.value(10),  # Adjust the x-position of the text box
    y=alt.value(10)   # Adjust the y-position of the text box
)
In [55]:
# Filter the data for the current year (2023)
qd_data_2023 = qd_monthly[qd_monthly['year'] == 2024]
qd_data_all = qd_monthly[qd_monthly['year'] != 2024]

# Create the boxplot for historical precipitation
qd_boxplot_TA = alt.Chart(qd_data_all).mark_boxplot(extent='min-max', size= 20, opacity=0.5).encode(
    alt.X('monthdate(date):O', title='Month (January-December)', axis=alt.Axis(format='%b')),  # Display month in "mmm" format
    alt.Y('rain:Q', title='Precipitation (mm)'),
).properties(width=500, height=250)

#Create a base chart for 'Mean_TA' values
qd_mean_ta_line = alt.Chart(qd_data_2023).mark_tick(size=20, thickness=3).encode(
    x=alt.X('monthdate(date):O', title=None, axis=alt.Axis(format='%b')),
    y=alt.Y('rain:Q', title='Precipitation (mm)'),
    tooltip=[alt.Tooltip('date', title="Date"), alt.Tooltip('rain', title="Precip")],
    color=alt.condition(
        alt.datum.percent > 1,
        alt.value("darkblue"),  # The positive color
        alt.value("darkorange")  # The negative color
    ),
).properties(
    width=500, height=250, title='Monthly Precipitation Totals'
)

# Create a chart for 'TA_Anomaly'
qd_anomaly_chart = alt.Chart(qd_data_2023).mark_bar().encode(
    x=alt.X('percent:Q', title=None, axis=alt.Axis(format='%'), scale=alt.Scale(domain=[0,1.8])),
    y=alt.Y('monthdate(date):O', title=None, axis=alt.Axis(format='%b')),
    color=alt.condition(
        alt.datum.percent > 1,
        alt.value("darkblue"),  # The positive color
        alt.value("darkorange")  # The negative color
    )
).properties(
    title = 'Percent of Normal',width=100, height=250
)

rule1 = alt.Chart(pd.DataFrame({
  'percent': ['1'],
  'color': ['gray']
})).mark_rule(opacity=0.8, strokeDash=[3,3]).encode(
  x='percent:Q',
  color=alt.Color('color:N', scale=None)
)
# Combine the two charts
text_percent = qd_anomaly_chart.mark_text(
    align='left',
    baseline='middle',
    size=10,
    dx=3  # Nudges text to right so it doesn't appear on top of the bar
).encode(
    text=alt.Text('percent:Q', format='.0%'),  # Format as percentage
)
anomaly_text = qd_anomaly_chart + text_percent

qd_combined_TA = qd_boxplot_TA + qd_mean_ta_line

combo_percent = rule1 + anomaly_text 

total = qd_combined_TA | combo_percent
total
Out[55]:
In [56]:
ranked = qd_pc[qd_pc['DOY'] == 365] 
ranked['rank'] = ranked['rain_cum_yr'].rank(ascending=False)
#ranked.head(22)

C:\Users\pmarshal\AppData\Local\Temp\ipykernel_5492\1804481392.py:2: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  ranked['rank'] = ranked['rain_cum_yr'].rank(ascending=False)

The following charts show cumulative annual precipitation. The first chart displays cumulative precipitation over the water year (Oct 1 - Sep 30), and the second chart shows Jan 1 to Dec 31. Total precipitation during the water year was very close to average; however, a very wet fall made the calendar year wetter than normal. 2024 was the sixth wettest year of the last 20 years at this weather station.

In [57]:
alt.data_transformers.disable_max_rows()

title_plot = alt.TitleParams(
   text='Lower Capilano Watershed - Cumulative Precipitation, Water Year',
   subtitle="Gray area is 5th to 95th percentile. Red dashed line is 20-year mean",
   anchor='middle',
   fontSize=14,
   fontWeight='bold')

pc_mean = alt.Chart(qd_pc[qd_pc['wtr_yr'] == 2022]).mark_line(color= "darkred", size= 1.5, strokeDash=[4, 2]).encode(
    alt.X('wtr_day_yr:T', title=None,
          axis=alt.Axis(format='%b-%d', tickCount=12, labelOverlap='parity')),
    alt.Y('mean:Q')
).properties(width=600)

area = alt.Chart(qd_pc[qd_pc['wtr_yr'] == 2022]).mark_area(color='lightgray', opacity=0.5).encode(
    alt.X('wtr_day_yr:T', title=None,
          axis=alt.Axis(format='%b-%d', tickCount=12, labelOverlap='parity')),
    alt.Y('percentile_95th:Q', title=''),
    alt.Y2('percentile_5th:Q')
).properties(width=600, height=300)

past_years = alt.Chart(last_few[last_few['wtr_yr'] != 2025], title=title_plot).mark_line(size= 2, opacity=0.8).encode(
    alt.X('wtr_day_yr:T', title=None,
          axis=alt.Axis(format='%b-%d', tickCount=12, labelOverlap='parity')),
    alt.Y('rain_cum_wyr:Q', title='Cumulative Precipitation (mm)'),
    color=alt.Color('wtr_yr:O', title='Water Year', scale=alt.Scale(scheme='set1')),
    size=alt.condition(
        alt.FieldEqualPredicate(field='wtr_yr', equal=2024),  # If year is 2025
        alt.value(3),  # Thicker line for 2025
        alt.value(2)),  # Default thickness for other years
    opacity=alt.condition(
        alt.FieldEqualPredicate(field='wtr_yr', equal=2024),  # If year is 2025
        alt.value(1),  # Thicker line for 2025
        alt.value(0.6)  # Default thickness for other years
    )
).properties(width=600)

total = area + pc_mean + past_years
total
Out[57]:
In [58]:
title3 = alt.TitleParams(
   text='Lower Capilano Watershed - Annual Cumulative Precipitation',
   subtitle="Gray area is 5th to 95th percentile. Red dashed line is 20-year mean",
   anchor='middle',
   fontSize=14,
   fontWeight='bold')
# Create a selection for the year
#year_selection = alt.selection_single(
#    name='Select',
#    fields=['year'],
#    bind=alt.binding_select(options=list(qd_pc['date'].dt.year.unique()), name='Select Year '),
#    init={'year': 2023}  # Initial selected year
#)

# Create a line chart using Altair
chart = alt.Chart(last_few[last_few['year'] != 2021]).mark_line().encode(
    alt.X('monthdate(date):T', title =None),
    alt.Y('rain_cum_yr:Q', title='Accumulated Precipitation (mm)'),
    color=alt.Color('year:N', title='Year', scale=alt.Scale(scheme='set1')),
    size=alt.condition(
        alt.FieldEqualPredicate(field='year', equal=2024),  # If year is 2025
        alt.value(3),  # Thicker line for 2025
        alt.value(2)),  # Default thickness for other years
    opacity=alt.condition(
        alt.FieldEqualPredicate(field='year', equal=2024),  # If year is 2025
        alt.value(1),  # Thicker line for 2025
        alt.value(0.6)  # Default thickness for other years
    )
).properties(
    width=600
)

area1 = alt.Chart(qd_pc[qd_pc['year'] == 2022]).mark_area(color='lightgray', opacity=0.5).encode(
    alt.X('monthdate(date):T', title=None),
    alt.Y('percent_95th:Q', title=''),
    alt.Y2('percent_5th:Q')
).properties(width=600, height=300)

# Create a line chart using Altair
mean_chart = alt.Chart(qd_pc, title=title3).mark_line(color='darkred', opacity=0.7, strokeDash=[4, 2], size= 1.5).encode(
    alt.X('monthdate(date):T', title =None),
    alt.Y('mean(rain_cum_yr):Q', title='Accumulated Precipitation (mm)')
).properties(
    width=600
)

area1 + chart + mean_chart
Out[58]:
In [59]:
mean_pc = qd_pc[(qd_pc['DOY'] == 365)]
mean_pc = mean_pc.iloc[:, 0:3]
mean_pc['average'] = mean_pc['rain_cum_yr'].mean()
mean_pc.tail()
percent = (2434/2864)*100
percent
Out[59]:
84.98603351955308
In [60]:
# Filter the data for the current year (2023)
qd_data_2023 = qd_monthly[qd_monthly['year'] == 2024]
qd_data_all = qd_monthly[qd_monthly['year'] != 2024]

qd_anomaly_chart = alt.Chart(qd_data_2023).mark_bar().encode(
    y=alt.X('percent:Q', title=None, axis=alt.Axis(format='%'), scale=alt.Scale(domain=[0,1.8])),
    x=alt.Y('monthdate(date):O', title=None, axis=alt.Axis(format='%b')),
    color=alt.condition(
        alt.datum.percent > 1,
        alt.value("darkblue"),  # The positive color
        alt.value("darkorange")  # The negative color
    )
).properties(width=500, height=250
)

rule1 = alt.Chart(pd.DataFrame({
  'percent': ['1'],
  'color': ['gray']
})).mark_rule(opacity=0.8, strokeDash=[3,3]).encode(
  x='percent:Q',
  color=alt.Color('color:N', scale=None)
)
# Combine the two charts
text_percent = qd_anomaly_chart.mark_text(
    align='left',
    baseline='middle',
    size=10,
    dx=-7,  # Nudges text to right so it doesn't appear on top of the bar
    dy=-7
).encode(
    text=alt.Text('percent:Q', format='.0%'),  # Format as percentage
)
anomaly_text = qd_anomaly_chart + text_percent

#anomaly_text

Atmospheric Rivers¶

The chart below highlights significant atmospheric river (AR) events in 2024. The October 18-20 event was a 1:50 to 1:100 year event in some locations. Heaviest rains fell at the southern end of the water supply areas, and in the municipalities closest to the North Shore Mountains. In total, more than 10% of the total annual precipitation fell during this one event. Total rainfall during this AR event equaled average rainfall amounts for the month of October.

In [61]:
title_1 = alt.TitleParams(
   text='Atmospheric Rivers in 2024',
   subtitle="Cumulative precipitation by event (Lower Capilano Watershed)",
   anchor='middle',
   fontSize=14,
   fontWeight='bold')
# Convert 'duration' and 'precip_cum' to float if they are not already
ar_df['duration'] = ar_df['duration'].astype(float)
ar_df['precip_cum'] = ar_df['precip_cum'].astype(float)

# Create the main line chart
ar_chart = alt.Chart(ar_df[ar_df['year'] == 2024], title=title_1).mark_line(opacity=0.8).encode(
    alt.X('duration:Q', title='Event Duration in Hours'),
    alt.Y('precip_cum:Q', title='Cumulative Precipitation (mm)'),
    alt.Color('monthdate(date):N', title='Date', sort='ascending'),
    tooltip=[alt.Tooltip('date', title="Date"), alt.Tooltip('precip_cum', title="Precip")]
).properties(width=500)


# Combine the main line chart, circle annotations, and text labels
ar_chart
Out[61]:
In [ ]: